! --------------------------------------------------------------------

subroutine parse_paramfile()
  use logs
  use myfiles
  use parameter_module
  use read_parameters
  implicit none

  call read_parameterfile(trim(adjustl(paramfile)))

#ifdef USE_CAM_FILE
  call get_parameter('CameraFile',params%CameraFile,required=.true.)
#endif
  call get_parameter('SnapDir',   params%SnapDir,   required=.true.)
  call get_parameter('SnapBase',  params%SnapBase,  required=.true.)
  call get_parameter('ObsDir',    params%ObsDir,    required=.true.)
  call get_parameter('ObsBase',   params%ObsBase,   required=.true.)
  call get_parameter('HsmlDir',   params%HsmlDir,   required=.true.)
  call get_parameter('HsmlBase',  params%HsmlBase,  required=.true.)
  call get_parameter('OutputDir', params%OutputDir, required=.true.)
  call get_parameter('VarFile',   params%VarFile,   required=.true.)
  call get_parameter('nx',        params%nx,        required=.true.)
#ifdef CUBE
  call get_parameter('ny',        params%ny,        required=.true.)
#endif
  call get_parameter('nz',        params%nz,        required=.true.)

  call finalize_parameters(abort_on_unused=.true.,abort_on_default=.false.)

end subroutine parse_paramfile

! --------------------------------------------------------------------

subroutine read_camerafile

#ifdef USE_CAM_FILE

  use logs
  use hdf5_wrapper
  use camera_module
  use parameter_module
#ifdef NONEVOLVING
  use arg_module, only : isnap,iframe
#else
  use arg_module, only : isnap
#endif
  implicit none
  integer*4 :: i
  integer*4 :: cam_handle,igood = 0

  call hdf5_open_file(cam_handle,trim(adjustl(params%camerafile)),readonly=.true.)
  call hdf5_read_attribute(cam_handle,'Flightpath/NumFrames',NumFramesCamFile)
  write(logunit,*)' -reading',NumFramesCamFile,'positions from Camerafile'
  call flush(logunit)

  call allocate_camera_arrays()
  call hdf5_read_data(cam_handle,'Flightpath/Frames'       ,cam_Frames)
  call hdf5_read_data(cam_handle,'Flightpath/X-coordinates',cam_x)
  call hdf5_read_data(cam_handle,'Flightpath/Y-coordinates',cam_y)
  call hdf5_read_data(cam_handle,'Flightpath/Z-coordinates',cam_z)
#ifdef VELOCITY
  call hdf5_read_data(cam_handle,'Flightpath/X-velocity',cam_vx)
  call hdf5_read_data(cam_handle,'Flightpath/Y-velocity',cam_vy)
  call hdf5_read_data(cam_handle,'Flightpath/Z-velocity',cam_vz)
#endif
  call hdf5_read_data(cam_handle,'Flightpath/Lx'  ,cam_Lx)
#ifdef CUBE
  call hdf5_read_data(cam_handle,'Flightpath/Ly'  ,cam_Ly)
#endif
  call hdf5_read_data(cam_handle,'Flightpath/Lz'  ,cam_Lz)
#ifdef USE_EIGENVECTORS
  call hdf5_read_data(cam_handle,'Flightpath/X-Prime',cam_xprime)
  call hdf5_read_data(cam_handle,'Flightpath/Y-Prime',cam_yprime)
  call hdf5_read_data(cam_handle,'Flightpath/Z-Prime',cam_zprime)
#else
  call hdf5_read_data(cam_handle,'Flightpath/Theta',cam_theta)
  call hdf5_read_data(cam_handle,'Flightpath/Phi'  ,cam_phi)
  call hdf5_read_data(cam_handle,'Flightpath/Psi'  ,cam_psi)
#endif
  call hdf5_close_file(cam_handle)

  ! Check camerafile contains necessary data
  do i=1,NumFramesCamFile
#ifdef NONEVOLVING
     if(cam_Frames(i) == iframe) then 
#else
     if(cam_Frames(i) == isnap) then 
#endif
        igood = 1
        exit      
     endif
  enddo
  if(igood /= 1) call abort('camerafile does not contain data for this snapshot')

  ! Set camera variables
  x0             = cam_x(i)
  y0             = cam_y(i)
  z0             = cam_z(i)
#ifdef VELOCITY
  vx0            = cam_vx(i)
  vy0            = cam_vy(i)
  vz0            = cam_vz(i)
#endif
  L(1)           = cam_Lx(i)
  L(2)           = cam_Ly(i)
  L(3)           = cam_Lz(i)
#ifdef USE_EIGENVECTORS
  xprime(:)      = cam_xprime(:,i)
  yprime(:)      = cam_yprime(:,i)
  zprime(:)      = cam_zprime(:,i)
#else
  theta          = cam_theta(i)
  phi            = cam_phi(i)
  psi            = cam_psi(i)
#endif

#endif

end subroutine read_camerafile

! --------------------------------------------------------------------

subroutine read_variables()
  use logs
  use parameter_module
  use smoothing_variables
  implicit none
  integer*4 :: ios, ivar, itype

  open(unit=10,file=trim(params%varfile),status='old',form='formatted',iostat=ios)
  if(ios /= 0) call abort(" -cannot open variables file: "//trim(params%varfile))
  do ivar=1,nvars_max
     read(10,*,iostat=ios)var_tmp(ivar)
     if(ios /= 0) exit
  enddo
  nvars = ivar-1

  ! Parse the PartType tag off and flag each variable's type
  do ivar=1,nvars 
     if(var_tmp(ivar)(1:1) == '/') var_tmp(ivar) = var_tmp(ivar)(2:len_trim(var_tmp(ivar)))

     ! Check we're not using velocity without VELOCITY set
#ifndef VELOCITY
     if(var_tmp(ivar)(11:18)=='Velocity') call abort("Cannot smooth velocity fields without engaging -DVELOCITY")
#endif

     selectcase(var_tmp(ivar)(1:9))
     case('PartType0') 
        if(var_tmp(ivar)(11:18)=='Velocity')then
           nvar_by_type(0)=nvar_by_type(0)+1
           var(0,nvar_by_type(0)) = trim(var_tmp(ivar))//'_X'
           nvar_by_type(0)=nvar_by_type(0)+1
           var(0,nvar_by_type(0)) = trim(var_tmp(ivar))//'_Y'
           nvar_by_type(0)=nvar_by_type(0)+1
           var(0,nvar_by_type(0)) = trim(var_tmp(ivar))//'_Z'
        else
           nvar_by_type(0)=nvar_by_type(0)+1
           var(0,nvar_by_type(0)) = trim(var_tmp(ivar))
        endif
     case('PartType1')
        if(var_tmp(ivar)(11:18)=='Velocity')then
           nvar_by_type(1)=nvar_by_type(1)+1
           var(1,nvar_by_type(1)) = trim(var_tmp(ivar))//'_X'
           nvar_by_type(1)=nvar_by_type(1)+1
           var(1,nvar_by_type(1)) = trim(var_tmp(ivar))//'_Y'
           nvar_by_type(1)=nvar_by_type(1)+1
           var(1,nvar_by_type(1)) = trim(var_tmp(ivar))//'_Z'
        else
           nvar_by_type(1)=nvar_by_type(1)+1
           var(1,nvar_by_type(1)) = trim(var_tmp(ivar))
        endif
     case('PartType2') 
        if(var_tmp(ivar)(11:18)=='Velocity')then
           nvar_by_type(2)=nvar_by_type(2)+1
           var(2,nvar_by_type(2)) = trim(var_tmp(ivar))//'_X'
           nvar_by_type(2)=nvar_by_type(2)+1
           var(2,nvar_by_type(2)) = trim(var_tmp(ivar))//'_Y'
           nvar_by_type(2)=nvar_by_type(2)+1
           var(2,nvar_by_type(2)) = trim(var_tmp(ivar))//'_Z'
        else
           nvar_by_type(2)=nvar_by_type(2)+1
           var(2,nvar_by_type(2)) = trim(var_tmp(ivar))
        endif
     case('PartType3') 
        if(var_tmp(ivar)(11:18)=='Velocity')then
           nvar_by_type(3)=nvar_by_type(3)+1
           var(3,nvar_by_type(3)) = trim(var_tmp(ivar))//'_X'
           nvar_by_type(3)=nvar_by_type(3)+1
           var(3,nvar_by_type(3)) = trim(var_tmp(ivar))//'_Y'
           nvar_by_type(3)=nvar_by_type(3)+1
           var(3,nvar_by_type(3)) = trim(var_tmp(ivar))//'_Z'
        else
           nvar_by_type(3)=nvar_by_type(3)+1
           var(3,nvar_by_type(3)) = trim(var_tmp(ivar))
        endif
     case('PartType4') 
        if(var_tmp(ivar)(11:18)=='Velocity')then
           nvar_by_type(4)=nvar_by_type(4)+1
           var(4,nvar_by_type(4)) = trim(var_tmp(ivar))//'_X'
           nvar_by_type(4)=nvar_by_type(4)+1
           var(4,nvar_by_type(4)) = trim(var_tmp(ivar))//'_Y'
           nvar_by_type(4)=nvar_by_type(4)+1
           var(4,nvar_by_type(4)) = trim(var_tmp(ivar))//'_Z'
        else
           nvar_by_type(4)=nvar_by_type(4)+1
           var(4,nvar_by_type(4)) = trim(var_tmp(ivar))
        endif
     case('PartType5') 
        if(var_tmp(ivar)(11:18)=='Velocity')then
           nvar_by_type(5)=nvar_by_type(5)+1
           var(5,nvar_by_type(5)) = trim(var_tmp(ivar))//'_X'
           nvar_by_type(5)=nvar_by_type(5)+1
           var(5,nvar_by_type(5)) = trim(var_tmp(ivar))//'_Y'
           nvar_by_type(5)=nvar_by_type(5)+1
           var(5,nvar_by_type(5)) = trim(var_tmp(ivar))//'_Z'
        else
           nvar_by_type(5)=nvar_by_type(5)+1
           var(5,nvar_by_type(5)) = trim(var_tmp(ivar))
        endif
     case default 
        call abort(' -invalid particle type: '//var_tmp(ivar)(1:9))
     endselect
  enddo

end subroutine read_variables

! --------------------------------------------------------------------

subroutine read_data(itype,get_hsml)
  use my_mpi
  use logs
  use parameter_module
  use arg_module, only : isnap
#ifdef VELOCITY
  use camera_module, only : vx0,vy0,vz0
#endif
  use myfiles
  use header
  use particles
  use grid
  use hdf5_wrapper
  use smoothing_variables
  use gadget3_bg_variables_module
  implicit none
  integer*4 :: i
  integer :: snap_handle,rank,dims(8),ierr,iPE
  integer ipart,ifile,ntotal,ivar,itype,strlen
  integer start2D(2), count2D(2), start1D(1), count1D(1), NPart_Master
  integer icount,ioffset
  integer(hid_t) :: file_id,error
  character(len=256) :: PartTypeName,varname
  character(len=256) :: snap_no_ext,LocalSnapshot
  character(len=256) :: obs_no_ext ,LocalObservables
  character(len=16)  :: csnap,cfile
  logical :: imulti,get_hsml, in_snap
#ifdef EOS_TEMP
  real*4, parameter :: T_EOS = EOS_TEMP
#endif
  integer*4        :: datasize
  character(len=32) :: datatype
  real*8, allocatable :: real8(:)

  write(logunit,*)'  -reading data from snapshot'
  call flush(logunit)

  ! Look for snapshot as single file; if not there look for the .0 etc extension 
  ! Turn off hdf5 error messages
  call h5eset_auto_f(0,error)  

  write(csnap,'(I3.3)')isnap 

#ifdef NO_SNAPDIRS
  snap_no_ext = trim(adjustl(params%snapdir))//'/'//trim(adjustl(params%snapbase))//'_'//trim(adjustl(csnap))
  obs_no_ext  = trim(adjustl(paramsobsdir))//'/'//trim(adjustl(params%obsbase))//'_'//trim(adjustl(csnap))
#else
  snap_no_ext = trim(adjustl(params%snapdir))//'snipshot_'//trim(adjustl(csnap))&
       //'/'//trim(adjustl(params%snapbase))//'_'//trim(adjustl(csnap))
  obs_no_ext = trim(adjustl(params%obsdir))//'observables_'//trim(adjustl(csnap))&
       //'/'//trim(adjustl(params%obsbase))//'_'//trim(adjustl(csnap))
#endif


  ! Test for presence of file - if not there, add on .0 to the filename
  LocalSnapshot = trim(adjustl(snap_no_ext))//'.hdf5'

  call h5fopen_f(LocalSnapshot,H5F_ACC_RDONLY_F,file_id,error)
  if(error.lt.0) then
     ! Can't open, try .0.hdf5 instead and try again
     LocalSnapshot = trim(adjustl(snap_no_ext))//'.0.hdf5'
     call h5fopen_f(LocalSnapshot,H5F_ACC_RDONLY_F,file_id,error)
     if(error < 0) then
        call abort('cannot find snapshot in single- or multi-file format: '//trim(adjustl(snap_no_ext)))
     else
        imulti = .TRUE.
        call h5fclose_f(file_id,error)
     endif
  else
     imulti = .FALSE.
     call h5fclose_f(file_id,error)
  endif

  ! Turn on hdf5 error messages again
  call h5eset_auto_f(1,error)

  ! Open this snapshot and read header to find how many files there are per snapshot
  call hdf5_open_file(snap_handle,LocalSnapshot, readonly=.true.)
  call read_header(snap_handle)
  call hdf5_close_file(snap_handle)

  ! Loop over the files and read in the particle numbers to fix the memory allocation
  NumPart_ByPE(:) = 0
  do ifile = 1,NumFilesPerSnapshot
     if(imulti)then
        write(cfile,*)ifile-1
        LocalSnapshot = trim(adjustl(snap_no_ext))//'.'//trim(adjustl(cfile))//'.hdf5'
     else
        LocalSnapshot = trim(adjustl(snap_no_ext))//'.hdf5'
     endif

     call hdf5_open_file(snap_handle,LocalSnapshot, readonly=.true.)
     call read_header(snap_handle)
     call hdf5_close_file(snap_handle)

     NumPart_Local = NumPart_ThisFile(itype)/NumPEs
     ntotal = NumPart_Local*(NumPEs-1)
     NPart_Master = NumPart_ThisFile(itype) - ntotal 
     if(MyPE == 0) NumPart_Local = NPart_Master
     NumPart_ByPE(MyPE) = NumPart_ByPE(MyPE) +  NumPart_Local
  enddo

  ! allocate particle data
  write(logunit,*)'   -PE',MyPE,' holds',NumPart_ByPE(MyPE),' particles'
  call flush(logunit)
  call allocate_read_particles(itype)

  ! Now we've allocated memory for data, read it in
  NumPart_ByPE(MyPE) = 0
  ioffset = 1
  write(logunit,*)'  -reading file NumFiles = ',NumFilesPerSnapshot
  call flush(logunit)
  do ifile = 1,NumFilesPerSnapshot

     if(imulti)then
        write(cfile,*)ifile-1
        LocalSnapshot    = trim(adjustl(snap_no_ext))//'.'//trim(adjustl(cfile))//'.hdf5'
        LocalObservables = trim(adjustl(obs_no_ext))//'.'//trim(adjustl(cfile))//'.hdf5'
     else
        LocalSnapshot    = trim(adjustl(snap_no_ext))//'.hdf5'
        LocalObservables = trim(adjustl(obs_no_ext))//'.hdf5'
     endif

     call hdf5_open_file(snap_handle,LocalSnapshot, readonly=.true.)
     call read_header(snap_handle)

     NumPart_Local = NumPart_ThisFile(itype)/NumPEs
     ntotal = NumPart_Local*(NumPEs-1)
     NPart_Master = NumPart_ThisFile(itype) - ntotal 
     if(MyPE == 0) NumPart_Local = NPart_Master
     NumPart_ByPE(MyPE) = NumPart_ByPE(MyPE) +  NumPart_Local

     ! read particle data
     start1D(1) = NPart_Master + (MyPE-1) * NumPart_Local + 1
     start2D(1) = 1
     start2D(2) = NPart_Master + (MyPE-1) * NumPart_Local + 1
     count1D(1) = NumPart_Local
     count2D(1) = 3
     count2D(2) = NumPart_Local

     if (NumPart_Local .gt. 0)then

        ! Read required variables first
        write(PartTypeName,'("PartType",I1)') itype

        varname = trim(PartTypeName)//'/Coordinates'
        call hdf5_read_data(snap_handle,varname,readCoordinates(1:3,ioffset:ioffset+NumPart_Local-1),&
             start=start2D,count=count2D)

#ifdef VELOCITY
        varname = trim(PartTypeName)//'/Velocity'
        call hdf5_read_data(snap_handle,varname,readVelocity(1:3,ioffset:ioffset+NumPart_Local-1),&
             start=start2D,count=count2D)
        
        ! Subtract offset velocity at source
        readVelocity(1,ioffset:ioffset+NumPart_Local-1) = readVelocity(1,ioffset:ioffset+NumPart_Local-1) - vx0
        readVelocity(2,ioffset:ioffset+NumPart_Local-1) = readVelocity(2,ioffset:ioffset+NumPart_Local-1) - vy0
        readVelocity(3,ioffset:ioffset+NumPart_Local-1) = readVelocity(3,ioffset:ioffset+NumPart_Local-1) - vz0
#endif

        ! Check header first - if non-zero, use that value to 
        ! fill array
        varname = trim(PartTypeName)//'/Mass'
        if(MassTable(itype+1) < 1.0d-20)then 
           call hdf5_read_data(snap_handle,varname,readMass(ioffset:ioffset+NumPart_Local-1),&
                start=start1D,count=count1D)
        else
           readMass(ioffset:ioffset+NumPart_Local-1) = real(MassTable(itype+1))
        endif

#ifdef FORCE_EOS
        if(itype == 0)then 
           varname = trim(PartTypeName)//'/OnEquationOfState'
           call hdf5_read_data(snap_handle,varname,readEOS(ioffset:ioffset+NumPart_Local-1),&
                start=start1D,count=count1D)
        endif
#endif

        do ivar=1,nvar_by_type(itype)

           if(trim(var(itype,ivar))=='PartType0/Weight') cycle
           if(trim(var(itype,ivar))=='PartType1/Weight') cycle
           if(trim(var(itype,ivar))=='PartType2/Weight') cycle
           if(trim(var(itype,ivar))=='PartType3/Weight') cycle
           if(trim(var(itype,ivar))=='PartType5/Weight') cycle
           if(trim(var(itype,ivar))=='PartType6/Weight') cycle
           if(trim(var(itype,ivar))=='PartType7/Weight') cycle

           ! Read the actual data: if velocity, we need to grab the values read earlier and manually enter
           ! them into three separate arrays
           selectcase(var(itype,ivar)(11:20))
#ifdef VELOCITY
           case('Velocity_X')
              readSmoothVar(ivar,ioffset:ioffset+NumPart_Local-1) = readVelocity(1,ioffset:ioffset+NumPart_Local-1)
           case('Velocity_Y')
              readSmoothVar(ivar,ioffset:ioffset+NumPart_Local-1) = readVelocity(2,ioffset:ioffset+NumPart_Local-1)
           case('Velocity_Z')
              readSmoothVar(ivar,ioffset:ioffset+NumPart_Local-1) = readVelocity(3,ioffset:ioffset+NumPart_Local-1)
#endif
           case default

              varname = trim(var(itype,ivar))

              ! Now need to check if the variable is on the list of standard Gadget3-BG variables.
              ! If not, we'll assume it's in the observables file rather than the snapshot
              in_snap = .FALSE.
              do i=1,n_bg_vars
                 if(varname(11:len_trim(varname)) == trim(bg_vars(i))) in_snap = .TRUE.
              enddo

              if(.not.(in_snap))then 
                 call hdf5_close_file(snap_handle)
                 call hdf5_open_file(snap_handle,LocalObservables, readonly=.true.)
              endif

              ! Check the precision of the smoothing variable
              call hdf5_get_type(snap_handle,varname,datatype,datasize)
              call flush(logunit)

              ! This is a complete kludge...
              if(datasize == 8)then 
                 allocate(real8(NumPart_Local))
                 call hdf5_read_data(snap_handle,varname,real8,&
                      start=start1D,count=count1D)                 
                 readSmoothVar(ivar,ioffset:ioffset+NumPart_Local-1) = real(real8(:)/1.0d30)
                 deallocate(real8)
              else
                 call hdf5_read_data(snap_handle,varname,readSmoothVar(ivar,ioffset:ioffset+NumPart_Local-1),&
                      start=start1D,count=count1D)
              endif

              ! Switch back to snapshot
              if(.not.(in_snap))then 
                 call hdf5_close_file(snap_handle)
                 call hdf5_open_file(snap_handle,LocalSnapshot, readonly=.true.)
              endif

           end select

#ifdef EOS_TEMP
           if(itype == 0)then 
              ! If smoothvar is temperature, set EOS particles to specified 
              if(varname == 'PartType0/Temperature')then 
                 do ipart=ioffset,ioffset+NumPart_Local-1
                    if(readEOS(ipart) > 0.5e0) readSmoothVar(ivar,ipart) = T_EOS
                 enddo
              endif
           endif

#endif
        enddo

        ! If get_hsml is set, read the smoothing lengths, otherwise read the ParticleIDs
        ! as we'll need to match them against the hsmlfile list
        if(get_hsml)then
            varname = trim(PartTypeName)//'/SmoothingLength'
           call hdf5_read_data(snap_handle,varname,readSmoothingLength(ioffset:ioffset+NumPart_Local-1),&
                start=start1D,count=count1D)
        else
           varname = trim(PartTypeName)//'/ParticleIDs'
           call hdf5_read_data(snap_handle,varname,readParticleIDs(ioffset:ioffset+NumPart_Local-1),&
                start=start1D,count=count1D)
        endif


        call hdf5_close_file(snap_handle)
        ioffset = ioffset + NumPart_Local

     endif
  enddo
  write(logunit,*)'   -done'
  call flush(logunit)

  ! Broadcast local particle number to each processor
  call mpi_barrier(mpi_comm_world,ierr) 
  do iPE=0,NumPEs-1
     call mpi_bcast(NumPart_ByPE(iPE),1,mpi_integer,iPE,mpi_comm_world,ierr)
  enddo

end subroutine read_data

! -------------------------------------------------------------------

subroutine read_data_catalogue(itype,get_hsml)
  use my_mpi
  use logs
  use parameter_module
  use arg_module, only : isnap
  use myfiles
  use header
  use particles
  use grid
  use hdf5_wrapper
  use smoothing_variables
  implicit none
  integer :: snap_handle,rank,dims(8),ierr,iPE
  integer ipart,ifile,ntotal,ivar,itype,strlen
  integer start2D(2), count2D(2), start1D(1), count1D(1), NPart_Master
  integer icount,ioffset
  integer(hid_t) :: file_id,error
  character(len=256) :: PartTypeName,varname
  character(len=256) :: infile_no_ext,LocalSnapshot
  character(len=16)  :: csnap,cfile
  logical :: imulti,get_hsml
#ifdef EOS_TEMP
    real*4, parameter :: T_EOS = EOS_TEMP
#endif


  write(logunit,*)'  -reading data from catalogue'
  call flush(logunit)

  ! Catalogues are limited to single file
  NumFilesPerSnapshot = 1

  ! Open the catalogue and read the numbers of particles
  open(unit=10,file=trim(adjustl(params%catalogue)),form='formatted',&
       status='old')
  read(10,*)NumPart_ThisFile
  read(10,*)BoxSize
  close(10)

  write(logunit,*)NumPart_ThisFile
  call flush(logunit)

  NumPart_ByPE(:) = 0
  NumPart_Local = NumPart_ThisFile(itype)/NumPEs
  ntotal = NumPart_Local*(NumPEs-1)
  NPart_Master = NumPart_ThisFile(itype) - ntotal 
  if(MyPE == 0) NumPart_Local = NPart_Master
  NumPart_ByPE(MyPE) = NumPart_ByPE(MyPE) +  NumPart_Local

  ! allocate particle data
  write(logunit,*)'   -PE',MyPE,' holds',NumPart_ByPE(MyPE),' particles'
  call flush(logunit)
  call allocate_read_particles(itype)

  ! Now we've allocated memory for data, read it in
  open(unit=10,file=trim(adjustl(params%catalogue)),form='formatted',&
       status='old')
  read(10,*) ! skip

  istart = NPart_Master + (MyPE-1) * NumPart_Local + 1
  iend   = istart + NumPart_Local - 1

  ioffset = 1
  if (NumPart_Local .gt. 0)then
     do ipart=1,NumPart_ThisFile(itype)
        if(ipart >= istart .and. ipart <= iend)then 
           read(10,*)readCoordinates(1,ioffset),readCoordinates(2,ioffset),readCoordinates(3,ioffset),&
                readMass(ioffset),readSmoothingLength(ioffset),readEOS(ioffset),readSmoothVar(1,ioffset)
           ! Advance counter
           ioffset = ioffset + 1
        else
           read(10,*) !Skip line
        endif
     enddo
  endif

  write(logunit,*)' -total(thispe) of Lx = ',sum(readSmoothVar(1,:)),'minmax',minval(readSmoothVar(1,:)),maxval(readSmoothVar(1,:))
  call flush(logunit)

  ! Broadcast local particle number to each processor
  call mpi_barrier(mpi_comm_world,ierr) 
  do iPE=0,NumPEs-1
     call mpi_bcast(NumPart_ByPE(iPE),1,mpi_integer,iPE,mpi_comm_world,ierr)
  enddo

end subroutine read_data_catalogue


! --------------------------------------------------------------------

subroutine read_hsmlfile(itype)
  use logs
  use my_mpi
  use parameter_module
  use arg_module, only : isnap
  use hdf5_wrapper
  use hsmlfile
  use seek_module
  use particles
  use smoothing_variables
  implicit none
  integer*4          :: itype,error,file_handle,ipart,ivar
  character(len=16)  :: csnap
  character(len=256) :: infile
  character(len=256) :: PartTypeName,varname

  write(logunit,*)'  -reading data from hsmlfile',itype
  call flush(logunit)

  write(csnap,'(I3.3)')isnap
  infile = trim(adjustl(params%hsmldir))//'/hsmldir_'//trim(adjustl(csnap))&
       //'/'//trim(adjustl(params%hsmlbase))//'_'//trim(adjustl(csnap))//'.hdf5'

  write(logunit,*)'  ',trim(infile)
  call flush(logunit)

  ! Allocate memory for ParticleIDs and SmoothingLength
  call allocate_hsmlfile()

  ! Open this snapshot and read header to find how many files there are per snapshot
  write(PartTypeName,'("PartType",I1)') itype
  call hdf5_open_file(file_handle,infile, readonly=.true.)
  varname = trim(PartTypeName)//'/ParticleIDs'
  call hdf5_read_data(file_handle,varname,hsml%ParticleIDs)
  varname = trim(PartTypeName)//'/SmoothingLength'
  call hdf5_read_data(file_handle,varname,hsml%SmoothingLength)
  call hdf5_close_file(file_handle)

  write(logunit,*)'   -done'
  call flush(logunit)

  write(logunit,*)'  -cross-referencing ParticleIDs from snapshot and hsmlfile'
  call flush(logunit)
  do ipart=1,NumPart_ByPE(MyPE)
     readSmoothingLength(ipart) = hsml%SmoothingLength(rb_seek(itype,readParticleIDs(ipart)))
  enddo

  write(logunit,*)'   -done'
  call flush(logunit)  

end subroutine read_hsmlfile

! --------------------------------------------------------------------

subroutine read_header(snap_handle)
  use header
  use hdf5_wrapper
  implicit none
  integer, intent(in) ::  snap_handle
  !
  call hdf5_read_attribute(snap_handle,'Header/NumPart_ThisFile',NumPart_ThisFile)
  call hdf5_read_attribute(snap_handle,'Header/NumPart_Total',NumPart_Total)
  call hdf5_read_attribute(snap_handle,'Header/NumFilesPerSnapshot',NumFilesPerSnapshot)
  call hdf5_read_attribute(snap_handle,'Header/Flag_Sfr',Flag_Sfr)
  call hdf5_read_attribute(snap_handle,'Header/Flag_Cooling',Flag_Cooling)
  call hdf5_read_attribute(snap_handle,'Header/Flag_StellarAge',Flag_StellarAge)
  call hdf5_read_attribute(snap_handle,'Header/Flag_Metals',Flag_Metals)
  call hdf5_read_attribute(snap_handle,'Header/Flag_Feedback',Flag_Feedback)
  call hdf5_read_attribute(snap_handle,'Header/MassTable',MassTable)

  call hdf5_read_attribute(snap_handle,'Header/Redshift',Redshift)
  call hdf5_read_attribute(snap_handle,'Header/BoxSize',BoxSize)
  call hdf5_read_attribute(snap_handle,'Header/Omega0',Omega0)
  call hdf5_read_attribute(snap_handle,'Header/OmegaLambda',OmegaLambda)
  call hdf5_read_attribute(snap_handle,'Header/HubbleParam',HubbleParam)

  ! Gadget3 features a different nomenclature
#ifdef GADGET2  
  call hdf5_read_attribute(snap_handle,'Header/Time',ExpansionFactor)
#else
  call hdf5_read_attribute(snap_handle,'Header/NumPart_Total_HighWord',NumPart_Total_HighWord)
  call hdf5_read_attribute(snap_handle,'Header/Time_GYR',Time_GYR)
  call hdf5_read_attribute(snap_handle,'Header/OmegaBaryon',OmegaBaryon)
  call hdf5_read_attribute(snap_handle,'Header/ExpansionFactor',ExpansionFactor)
#endif


end subroutine read_header

! --------------------------------------------------------------------

subroutine write_header(snap_handle)
  use header
  use hdf5_wrapper
  implicit none
  integer, intent(in):: snap_handle

  call hdf5_create_group(snap_handle,"/Header")
  !
  call hdf5_write_attribute(snap_handle,'Header/NumPart_ThisFile',NumPart_ThisFile)
  call hdf5_write_attribute(snap_handle,'Header/NumPart_Total',NumPart_Total)
  call hdf5_write_attribute(snap_handle,'Header/NumPart_Total_HighWord',NumPart_Total_HighWord)
  call hdf5_write_attribute(snap_handle,'Header/NumFilesPerSnapshot',NumFilesPerSnapshot)
  call hdf5_write_attribute(snap_handle,'Header/Flag_Sfr',Flag_Sfr)
  call hdf5_write_attribute(snap_handle,'Header/Flag_Cooling',Flag_Cooling)
  call hdf5_write_attribute(snap_handle,'Header/Flag_StellarAge',Flag_StellarAge)
  call hdf5_write_attribute(snap_handle,'Header/Flag_Metals',Flag_Metals)
  call hdf5_write_attribute(snap_handle,'Header/Flag_Feedback',Flag_Feedback)
  call hdf5_write_attribute(snap_handle,'Header/MassTable',MassTable)
  call hdf5_write_attribute(snap_handle,'Header/Redshift',Redshift)
  call hdf5_write_attribute(snap_handle,'Header/BoxSize',BoxSize)
  call hdf5_write_attribute(snap_handle,'Header/Omega0',Omega0)
  call hdf5_write_attribute(snap_handle,'Header/OmegaLambda',OmegaLambda)
  call hdf5_write_attribute(snap_handle,'Header/HubbleParam',HubbleParam)
#ifdef GADGET2  
  call hdf5_write_attribute(snap_handle,'Header/Time',ExpansionFactor)
#else
  call hdf5_write_attribute(snap_handle,'Header/ExpansionFactor',ExpansionFactor)
  call hdf5_write_attribute(snap_handle,'Header/Time_GYR',Time_GYR)
  call hdf5_write_attribute(snap_handle,'Header/OmegaBaryon',OmegaBaryon)
#endif

end subroutine write_header
